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Prediction of the backflow and 
recovery regions in the backward facing 
step at various Reynolds numbers 

By V. Michelassi 1 5 P. A. Durbin 2 AND N. N. Mansour 3 


A four-equation model of turbulence is applied to the numerical simulation of flows 
with massive separation induced by a sudden expansion. The model constants are 
a function of the flow parameters, and two different formulations for these functions 
are tested. The results are compared with experimental data for a high Reynolds- 
number case and with experimental and DNS data for a low Reynolds-number 
case. The computations prove that the recovery region downstream of the massive 
separation is properly modeled only for the high Re case. The problems in this case 
stem from the gradient diffusion hypothesis, which underestimates the turbulent- 
diffusion. 


1. Introduction 

The Reynolds Averaged Navier Stokes equations (RANS) equations need a tur- 
bulence model for computation of Reynolds stresses that stem from averaging the 
non-linear convective terms. A large family of turbulence models exists in the liter- 
ature. The models range from simple algebraic expressions for the eddy viscosity to 
more elaborate formulations which introduce a separate transport equation for each 
component of the Reynolds Stress tensor. Eddy viscosity models such as the k - e 
model still represent a good compromise between accuracy and computational effi- 
ciency and will be the subject of this investigation. Moreover, the results of a recent 
workshop (Rodi et al , 1995) showed that, even though full Reynolds stress models 
bring more physics into the model, the large increment in the computational effort 
associated with these models is not always followed by a proportional improvement 
in the quality of the predictions. 

Two-equation models of turbulence have been recently tuned with the aid of 
Direct Numerical Simulation (DNS) data (see e.g . Michelassi and Shih, 1991, Rodi 
et a/., 1993). This tuning was mostly done to allow modeling of the near wall 
region and to reproduce the profiles of the turbulent kinetic energy, k, and of the 
dissipation rate e in this critical flow region. The tuning was done by using fully 
developed or turbulent boundary layer flows (Rodi and Mansour, 1990). Most of the 
so called “low Reynolds number modifications” ( LR ) to the two-equation models of 
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turbulence were able to improve the model capability in the flow layer close to the 
wall. Nevertheless, little or no change at all was found in the core region of the flow 
since most of the modifications were designed to vanish away from solid boundaries 
(Zhu and Shih, 1993). 

The LR models, which allow the integration of the equations in the near wall 
region, can successfully model a wide range of flows, but often do not allow flows 
with strong adverse pressure gradients and/or separation to be computed accurately. 
This seems to be a general problem associated with the two-equation formulation 
irrespective of the treatment of the near-wall region (Michelassi, 1993). In the 
backward facing step flow, both an adverse pressure gradient and flow separation 
are to be modeled, which makes this test case particularly challenging. 

Durbin (1995) computed the backward facing step flow at different Reynolds 
numbers. His computations proved that downstream of the reattachment point the 
computed velocity profiles tend too slowly to a boundary layer profile for the high 
Reynolds number case, but not for the low Reynolds number case. A similar failure 
was encountered by Rodi (1991) with a two-layer model of turbulence. Again, 
the velocity profiles in the recovery region tend too slowly to a developed profile. 
Durbin and Rodi use forms of the two-equation k — e model which, while based on 
the Boussinesq assumption, have very little in common with the treatment. This 
indicates that the problems are stemming from the k — e frame and not from the 
wall treatment. 

This phenomenon is also of great importance in practical flows with engineering 
relevance such as the flow in turbomachines. In fact, immediately downstream of 
the trailing edge of a turbine or a compressor blade, two counterrotating vortices 
interact with the wake in a very similar manner to that found for the backward 
facing step. The modeling of the wake downstream of the two vortices is of pri- 
mary importance in turbomachinery flows because of its impact on the stator-rotor 
interaction. In this case, the computed wake decay, which is similar to the flow 
recovery region in the backward facing step, seems to be too slow compared to the 
measurements as indicated by a number of computations for subsonic and transonic 
turbines (Michelassi et al 1995). These results were shown to be true regardless of 
the assumption of a fully turbulent or transitional boundary layer along the blade 
profile. In the turbomachinery flow case, it is not clear if the discrepancies are due 
to the inherently unsteady nature of the experimental flow field, or to deficiencies 
in the model as in the backward facing step where the steadiness of the flow is not 
an issue. 

Although the recovery region problem with computing the backstep has been 
often pointed out, very little has been done so far to identify the causes of the slow 
recovery downstream of the reattachment point. Two-equation models are known to 
have theoretical limitations which stem mainly from the eddy viscosity assumption. 
Still, the ability of these simple turbulence models to mimic a flow with massive 
separation and the wake decay needs to be improved. 

With this in mind three different backward facing step data sets are used to 
compare with the computations and to identify the reasons for the discrepancies 
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between computations and measurements in the recovery region. 


2. The turbulence model 

The turbulence model uses the standard k — e equations: 

d t k + U • Vfc = P k -£ + [( 1 / + — )Vk], (1) 

<?k 

0t£ + U - Ve = Cf ' Pk ~ — + [(v + — )Vf]. (2) 

1 Vt 

The model constant C ( i is computed as: 
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in which cl is the minimum distance from the wall, and L is the turbulence length 
scale. On no-slip boundaries, y — ► 0, 


k = 0, 


2v 


k 


Two additional equations are solved. The first transport equation determines the 
velocity fluctuation normal to the wall, v 2 . The v 2 transport equation is 


d t v 2 + U ■ Vv 2 = kf -v 2 j + V- l(t' + i/,)Vu2], 


(4) 


where kf represents redistribution of turbulence energy from the streamwise com- 
ponent. Non-locality is represented by solving an elliptic relaxation equation for 
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in which 
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The Boussinesq approximation is used for the stress-strain relation: 


_ «i«> 2 _ 

aij - ~r ~ 3 ~~k ,j ' 


where the eddy viscosity is given by 
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The constants of the model are: 


= 0.19, a k = 1, a ( = 1.3, C (1 = 1.55, C (2 = 1.9 
C l = 1.4, C 2 = 0.3, C L = 0.3, C„ = 70. 


The boundary conditions are 


( 7 ) 


v 2 = 0 , /( 0 ) 


20u 2 v 2 

e(0)y 4 ' 


on no-slip walls. 

The original model formulation was modified by Durbin and Laurence (1996) 
in the expressions for the length and time scales, L and T, and the definition of 
the model constant C ( \. The length and time scales are now computed to allow a 
smoother switch from the core-flow values to the near- wall values as follows: 

L 2 = C 2 (k 3 /e 2 + Cy' 2 /e'/ 2 ) , T 2 = k 2 /e 2 + C 2 v/e (8) 

The selected values of the constants are C p = 0.2, C ^ — 70, and Ct = 6. 

In Eq. (3) the scaling of C t \ in the near wall region is done by using the wall 
distance y. The definition of the wall distance can be problematic in complex flows 
so that Durbin and Laurence (1996) replaced Eq. (3) with another expression based 
on v 2 which is suited to feel the proximity of the wall: 

C ( i = 1.44 (l + 1/30 (fc/^) 1/2 ) , (9) 

This expression, like the one in Eq. (3), is supposed to increase the production of 
dissipation in the near wall region, where v 2 goes to zero faster than k. Both the 
original formulation, hereafter referred as form (1) of the model, and the modified 
formulation, hereafter referred as form (2), have been applied with the same inlet 
and boundary conditions. 

3. The data sets and the computations 

The turbulence model with the two different forms described in the previous 
section was applied to the computation of three different backward-facing step ge- 
ometries and different Reynolds numbers. 

The first experimental data set considered here is that of Driver and Seegmiller 
( DS ) (1985) which allowed testing the model in a high Reynolds-number configu- 
ration with a Reynolds number based on the step height of 37, 500. Measurements 
were taken by using laser velocimetry and include mean and instantaneous quanti- 
ties and triple correlations. 

The low Reynolds-number case refers to the measurements by Kasagi and Mat- 
sunaga (KM) (1995). In this case the flow Reynolds number, based again on the 
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Figure 1 . DS Velocity profiles, o experiments, model version (1), 

model version (2). 

step height, is 5540. Measurements were taken by using a particle image velocime- 
t.ry method (PIV) which allowed measuring instantaneous and average quantities. 
The measured profiles were also carefully tested to verify mass conservation. A 
similar Reynolds number (Re = 5100) was achieved by Le and Moin (LM) (1994) 
which produced a DNS data set for the backward facing step geometry. The large 
amount of information on the flow field makes this DNS data set very valuable for 
testing and developing two-equation models of turbulence. 

The investigation is carried out on three different data sets to test the model 
under different Reynolds number conditions. At the present stage of research it is 
still impossible to perform the DNS of a backward facing step at high Reynolds 
number, so the use of an experimental data set was compulsory. The two data sets 
for the low Reynolds number case were selected to verify that model testing done 
by using a classical experimental data set could be extended to the DNS data for 
such a flow field. 

The computational grids for the three test cases have 120 x 120 grid nodes clus- 
tered near solid walls. The inlet section profiles have been carefully specified as 
follows. For the DS case the inlet profiles have been computed by a boundary layer 
code until the momentum thickness of 5000 was reached (Durbin, 1995). These pro- 
files were then imposed at the inlet section of the computational domain. For the 
KM case the inlet profiles were those of a fully developed channel flow, as indicated 
by Kasagi and Matsunaga (1995) in their discussion of the flow nature upstream of 
the separation point. For the test case proposed by Le and Moin, the inlet profiles 
were those computed by the DNS at the section upstream of the separation point 
corresponding to the inlet section of the present computational grid. No other grids 
were used for the calculations since the 120 x 120 grid was already found adequate 
for this kind of computation by Durbin (1995). 

The first set of computations refer to the DS case. Fig. 1 compares the measured 
profiles with those computed by using the two versions of the model. In all the 
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FIGURE 2. DS Turbulent kinetic energy profiles. Symbols as in Fig. 1. 



FIGURE 3. DS Turbulent shear stress profiles. Symbols as in Fig. 1. 

following plots the ordinate y = 1 corresponds to the step corner. The reattachment 
point is not affected by the change in the model, but the different functions adopted 
for the computation of the length scale Z, the time scale T, and the coefficient of the 
production rate of dissipation C t \ show some effect in the backflow region. Here 
version (2) of the model moves the computed profiles closer to experiments. A 
sensitivity analysis made by changing the coefficients in Eqs. (8,9) proved that the 
model is sensitive to the value of C p , which was set equal to 0.2. The model can 
be seen to predict velocity profiles which are steeper than the measured ones in the 
backflow region. Moreover, in the recovery region the computations lag behind the 
experimental boundary layer profile. The agreement is indeed quite good in terms 
of turbulent kinetic energy (see Fig. 2) and turbulent shear stress (see Fig. 3). 
Apparently, the models succeed in reproducing the correct level of turbulent kinetic 
energy and shear stress with the only exception of a narrow region deep inside of 
the backflow, where the maximum of turbulent kinetic energy and turbulent shear 
stress are not correctly predicted and somewhat misplaced. 
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FIGURE 4. KM Velocity profiles. experiments, model version (1), 

model version (2). 



FIGURE 5. LM Velocity profiles. Symbols as in Fig. 4. 

When moving to the KM and LM test cases, a more careful analysis is possible 
due to the large number of measurements. Figure 4 compares the measured and 
computed velocity profiles in several stations starting from the separation point for 
the KM test case. The agreement is again quite good, and apparently the two 
versions of the model give almost identical results in this case. The recovery region 
is well predicted here. Again, the backflow region shows the steep velocity profiles 
predicted in the high Reynolds number case, while the experiments show a profile 
which seems to indicate quite a low turbulence level. Figure 5 shows the same veloc- 
ity profiles for the LM test case. In this last computation the recirculation bubble 
length was underestimated by approximately 4%. The backflow region length was 
computed in almost perfect agreement with the experiments for the I\M case. The 
plots also show that the differences in the computation of the length and time scales 
in the two versions of the code bring very little change to the computed profiles, 
which are almost collapsing on each other, in the low Re case. 




80 


V. Michelassi et al 



FIGURE 6. LM Turbulent kinetic energy profiles. Symbols as in Fig. 4. 



FIGURE 7. LM Turbulent shear stress profiles. Symbols as in Fig. 4. 

Figures 4 and 5 show that there is very little difference between the LM and 
KM data sets. Since the information given by the KM and LM cases do not show 
significant differences, only the latter will be described in detail in what follows. 

Figure 6 compares the measured and computed turbulent kinetic energy profiles 
at several stations starting from the separation point. The agreement between com- 
putations is generally satisfactory, even though the models overpredict the turbulent 
kinetic energy in the backflow region. Of the two, version (2) of the model seems to 
reduce the overprediction. This was also found in the high Reynolds number case. 
This overprediction spreads in the shear layer as the flow proceeds downstream. 

The overprediction of k seems to have an effect in terms of turbulent shear stress 
also, as shown in Fig. 7. Here the turbulent shear stress is overestimated by both 
the formulations in the backflow region and underestimated in the recovery region. 
The change from overestimation to underestimation takes place gradually across 
the reattachment point and the fit between DNS and computations improves only 
far downstream. 



Backward facing step flow 


81 



FlGliRK 8. LM Dissipation rate profiles. Symbols as in Fig. 4. 

The LM data set also includes the dissipation rate. Figure 8 shows that the 
computed dissipation rate level is larger than that given by the DNS in the backflow 
region. The e levels are well predicted in the recovery region. 

The skin friction coefficients in Figs. 9 and 10 show that version (2) of the model 
tends to reduce the recirculation bubble length in the low- Re number case (a similar 
trend was also found for the KM test), whereas the same model seems to increase 
the backflow region length in the high-Re case. In the KM case, also, a reduction 
of the recirculation bubble length was observed. 

4. Discussion and conclusions 

The brief description of the computations done in the previous section evidences 
how the computed overall flow pattern agrees with the high-i?e and low -Re cases, 
although some discrepancies between the computations and the measurements (and 
DNS) arise in terms of turbulence quantities. 

In the recovery region, as already pointed out by several authors, ( e.g . Durbin, 
1995), the computations recover to a boundary layer profile much more slowly than 
experiments would indicate at high Reynolds numbers. This disagreement fades 
away for smaller Reynolds numbers, as those typical of the DNS. In the backflow 
region the computed profiles seem too steep, which would indicate too large a tur- 
bulence level. 

Version (2) of the model was found to work slightly better than the original 
version of the model in the backflow region. This can be attributed to the different 
choice of the length scale formula. In version (1) the model chooses between two 
different values of the length scale, whereas in version (2) the expression for the 
length scale allows a smooth switch from the two values. Observe that the same 
smooth switch is guaranteed for the computation of the time scale. This seems to 
play a significant role in the improvement of the results where, due to the small 
local Reynolds number, the expressions for the length and time scale are switching 
between the two values. In the recovery region the local Reynolds number is larger 



82 


V. Michelassi et al 


and the beneficial effects of the smooth transition between the two values of the 
time and length scale formulas disappear. 

In terms of turbulent shear stress, the backflow region again shows some slight 
inaccuracies for both the high -Re and low-i?e cases. This fits with the shape of the 
computed velocity profile, which indicates that the mean velocity gradient and the 
turbulence levels are too high. From the DNS data set it is possible to compute 
a turbulent viscosity fi t via the definition of the turbulent shear stress given in 
the Boussinesq assumption. This sort of computation does not guarantee that the 
turbulent viscosity is positive, since there is no guarantee that the mean shear 
and the turbulent shear stress always have opposite sign: in fact Fig. 11 shows that 
turbulent viscosity computed from the turbulent shear stress by DNS gives negative 
values. 

The turbulent viscosity is small deep inside the backflow region and grows toward 
the reattachment point. The two versions of the model are found to overestimate 
the turbulent viscosity in the backflow region. There is very little difference be- 
tween the computations all through the computational domain. Observe that a 
large turbulent viscosity would imply a large momentum diffusion, which should 
decrease the recirculation bubble length. Surprisingly, this is not the case in the 
computations: the overestimation in (it is followed by an excellent agreement be- 
tween the computed and measured reattachment point. The figure also shows that 
the disagreement between the computations and the DNS fades away downstream 
of the reattachment point. But the same figure also shows that in the recovery 
region the turbulent viscosity is underestimated. The smaller momentum diffusion 
in the computation could partially explain why the computed velocity profiles tend 
to the boundary layer profiles too slowly. The discrepancies between DNS and com- 
putations are mainly in the backflow region and the shear layer, since above the 
latter the computations seem to follow the DNS profiles quite well. 

The DNS data set also contains all the terms of the transport equation for the 
turbulent kinetic energy. With these data it was possible to evaluate the accuracy 
of each term of the modeled transport equation for k. A full comparison of all the 
terms (i,e. convection = C\ t, viscous diffusion= V d, turbulent diffusion^ Td, pres- 
sure diffusion= Pd , production = P*, and dissipation = e) showed that the viscous 
diffusion Vd has nearly no effect. The computed convection of fc, CV, is in very good 
agreement with both the measurements by KM and the DNS by LM . The dissipa- 
tion rate, although not in perfect agreement with the data in the backflow region, 
closely resembles the DNS profiles in the shear layer. So, the terms which need a 
further check, and that are not often compared with the experiments for models 
based on the eddy viscosity, are the production rate P \ and both the turbulent and 
pressure diffusion terms Td and Pd respectively. 

Figure 12 compares the DNS production rate versus the profiles obtained by using 
the two different versions of the model. The agreement between computations and 
the DNS profiles is good. Observe that the peak in the production rate, which is 
probably caused by the very high mean shear downstream of the separation point, 
is well captured. The production rate is somewhat overpredicted in the backflow 
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FIGURE 9. DS Skill friction coefficient, o experiments, model version (1), 

model version (2). 



FIGURE 10. LM Skill friction coefficient. Symbols as in Fig. 9. 



FIGURE 11. LM Turbulent viscosity profiles. Symbols as in Fig. 9. 
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Pk + x 

FIGURE 12. LM Production rate profiles. experiments, model version 

(1), model version (2). 



Figure 13. LM Td + Pd profiles. Symbols as in Fig. 12. 

region, but this overprediction seems to fade away as the reattachment point is 
reached. The same agreement was found in the high -Re case. 

Before comparing the turbulent and pressure diffusion terms, one should recall 
that the gradient diffusion hypothesis, done in the k — e model, does not distinguish 
between Pd and Td , which are just lumped together. Still, it is possible to compare 
the sum oi Td and Pd from the DNS calculations with the computed turbulent 
diffusion of turbulent kinetic energy, which should be the sum of the two. Observe 
that the comparison is done for the diffusive terms (second order derivative of k for 
the k — e model and first order derivative of Td and Pd for the DNS data). Figure 
13 compares the computed diffusion of k with the sum of the turbulent diffusion and 
pressure diffusion from the DNS. The agreement between computations and DNS is 
quite good. The up-down shape of the profile from the DNS is closely reproduced by 
the calculations. The agreement remains good in the entire computational domain 
and does not deteriorate when making the same comparison for the KM data set. 
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FIGURE 14. DS Td profiles. Symbols as in Fig. 12. 

When making the same comparison by using the DS data base at a higher 
Reynolds number, some problems arise due to the scatter of the measured data. 
Figure 14 compares the turbulent diffusion of turbulent kinetic energy for the DS 
case. Although there are not as many data as in the DNS case, the figure clearly 
suggests that the turbulent diffusion is largely underestimated in the shear layer 
from the separation point till far downstream. The underestimation is quite severe 
and clearly limited to the flow region where the mean shear is high. However, the 
experimental data are probably not accurate enough to differentiate, as in Fig. 14. 

Figures 13 and 14 indicate that as long as the Reynolds number is small, the 
gradient diffusion hypothesis gives the correct estimate of the turbulent plus the 
pressure diffusion, especially in the high shear layer. The two figures also show that 
the same closure hypothesis fails when the Reynolds number is large. Apparently at 
large Re there is a large scatter of turbulence time and length scales. This scatter 
is probably not modeled when using a linear eddy viscosity model. The scatter is 
reduced at smaller Reynolds number, and the turbulence model then agrees much 
better with the experiments and DNS. 

In conclusion, the computations show that the slow recovery downstream of the 
reattachment point occurs only in high Reynolds number flows and is probably 
caused by the gradient diffusion hypothesis, which is not able to model the large 
turbulent diffusion typical of the high shear layer. In the backflow region the com- 
putations and the comparison with the experiments and the DNS do not allow 
identification of any specific deficiency of the model. Still, the plots indicate that 
in the backflow region the models predict too high a turbulence level and too much 
velocity gradient, which are interrelated deficiencies. 
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